rm(list=ls()) 
### Import datasets
all1ay <- read.csv("/Users/shannoncarcelli/Box Sync/UCSD/Dissertation Chapters/Crossnational/JOP Final Materials/all1ay.csv")
all1 <- read.csv("/Users/shannoncarcelli/Box Sync/UCSD/Dissertation Chapters/Crossnational/JOP Final Materials/all1.csv")
num <- aggregate(year ~ countryname, data = all1, FUN = function(x) length(unique(x)))
all <- subset(all1, countryname %in% subset(num$countryname, num$year == max(num$year)) == T)
allay <- subset(all1ay, countryname %in% subset(num$countryname, num$year == max(num$year)) == T)
### Delete loan-granting and regional agencies from agency-year datasets
all2ay <- subset(all1ay, type != "DFI" & type != "ExportBank" & type != "Regional")
aylimit <- subset(allay, type != "DFI" & type != "ExportBank" & type != "Regional")


#### REPLICATION CODE FOR TABLES AND FIGURES ###
### TABLES
## Table 1: State-Level Results
install.packages("multiwayvcov", "sandwich", "lmtest", "stargazer")
library(sandwich)
packageVersion("sandwich") ##Version must be ‘3.0.2’ or higher
library(lmtest)
library(stargazer)
## Plurality as DV
# Run regressions
f1p <- lm(frag ~ pluralty, data = all1)
f2p <- lm(frag ~ pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f3p <- lm(frag ~ factor(year) + pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f4p <- lm(frag ~ countryname + pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f5p <- lm(frag ~ factor(year) + countryname + pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
# Calculate (clustered) standard errors
vcov1fp <- vcovCL(f1p, cluster = all1$countryname)
vcov2fp <- vcovCL(f2p, cluster = all1$countryname)
vcov3fp <- vcovCL(f3p, cluster = all1$countryname)
vcov4fp <- vcovCL(f4p, cluster = all1$countryname)
vcov5fp <- vcovCL(f5p, cluster = all1$countryname)
# Re-run with new standard errors
f1cp <- coeftest(f1p, vcov = vcov1fp)
f2cp <- coeftest(f2p, vcov = vcov2fp)
f3cp <- coeftest(f3p, vcov = vcov3fp)
f4cp <- coeftest(f4p, vcov = vcov4fp)
f5cp <- coeftest(f5p, vcov = vcov5fp)
# Results in LaTex format
stargazer(f1cp, f2cp, f3cp, f4cp, f5cp, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Plurality", "Veto Players", "Lib. Econ.", "Gov. Vote-Share", "Opp. Vote-Share", "HHI Government", "Fractionalization", "Leg. Election Year", "Exec. Election Year", "Unified Gov.", "Polarization", "Conservative Gov.", "GDP (trillion USD)", "GDP/cap (100,000 USD)", "Aid Budget (100 billion USD)"), 
          dep.var.labels = "Fragmentation (1 - HHI)",
          add.lines = list(c("Year FE", "N", "N", "Y", "N", "Y"),
                           c("Donor FE", "N", "N", "N", "Y", "Y"),
                           c("Observations", paste(length(f1p$fitted.values)), paste(length(f2p$fitted.values)), paste(length(f3p$fitted.values)), paste(length(f4p$fitted.values)), paste(length(f5p$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(f1p)[9]), digits = 3)), paste(round(as.numeric(summary(f2p)[9]), digits = 3)), paste(round(as.numeric(summary(f3p)[9]), digits = 3)), paste(round(as.numeric(summary(f4p)[9]), digits = 3)), paste(round(as.numeric(summary(f5p)[9]), digits = 3)))))

## Table 2: Bureaucracy-Level Results
# Rescale "percent tied" DV for easier reading
all2ay$pct.rescale <- all2ay$percent2*100
all2ay$don_factor <- as.factor(all2ay$donor)
# Run regressions
tt1 <- lm(pct.rescale ~ main, data = all2ay)
tt2 <- lm(pct.rescale ~ main + pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all2ay)
tt4 <- lm(pct.rescale ~ factor(year) + main + pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all2ay)
tt5 <- lm(pct.rescale ~ donor + main + pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all2ay)
tt6 <- lm(pct.rescale ~ factor(year) + donor + main + pluralty + checks + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all2ay)
# Calculate (clustered) standard errors
vcov1 <- vcovCL(tt1, cluster = all2ay$don_factor)
vcov2 <- vcovCL(tt2, cluster = all2ay$don_factor)
vcov4 <- vcovCL(tt4, cluster = all2ay$don_factor)
vcov5 <- vcovCL(tt5, cluster = all2ay$don_factor)
vcov6 <- vcovCL(tt6, cluster = all2ay$don_factor)
# Re-run with new standard errors
tie1ayc <- coeftest(tt1, vcov = vcov1)
tie2ayc <- coeftest(tt2, vcov = vcov2)
tie4ayc <- coeftest(tt4, vcov = vcov4)
tie5ayc <- coeftest(tt5, vcov = vcov5)
tie6ayc <- coeftest(tt6, vcov = vcov6)
# Results in LaTex format
stargazer(tie1ayc, tie2ayc, tie4ayc, tie5ayc, tie6ayc, no.space = T, digits = 2, omit = c("year", "donor", "dy", "type"), 
          covariate.labels = c("Primary Agency", "Plurality", "Veto Players", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          add.lines = list(c("Year FE", "N", "N", "Y", "N", "Y"),
                           c("Donor FE", "N", "N", "N", "Y", "Y"),
                           c("Observations", paste(length(tt1$fitted.values)), paste(length(tt2$fitted.values)), paste(length(tt4$fitted.values)), paste(length(tt5$fitted.values)), paste(length(tt6$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(tt1)[9]), digits = 3)), paste(round(as.numeric(summary(tt2)[9]), digits = 3)), paste(round(as.numeric(summary(tt4)[9]), digits = 3)), paste(round(as.numeric(summary(tt5)[9]), digits = 3)), paste(round(as.numeric(summary(tt6)[9]), digits = 3)))))

### FIGURES
## Figure 1: Mean Bureaucratic Fragmentation of Aid, European Donors
frag <- aggregate(frag ~ countryname, data = all1, FUN = function(x) mean(x, na.rm = T))
frag$countryname <- ifelse(frag$countryname == "FRG/Germany", "Germany", paste(frag$countryname))
frag1<- subset(frag, frag > 0 & (countrycode(countryname, "country.name", "region23") == "Western Europe" | countrycode(countryname, "country.name", "region23") == "Northern Europe"))
#frag1<- subset(frag, frag > 0 & (countrycode(countryname, "country.name", "region23") == "Western Europe" | countrycode(countryname, "country.name", "region23") == "Northern Europe"| countrycode(countryname, "country.name", "region23") == "Northern America"))
par(mar=c(4,6,2,1)+.4)
frag2 <- frag1[order(frag1$frag),]
barplot(frag2$frag, names.arg = frag2$countryname, horiz = T, cex.names = 1, las = 2, main = "Mean Bureaucratic Fragmentation of Aid, \nEuropean Donors", xlab = "1 - Herfindahl-Hirshman Index", col = "skyblue4")

## Figure 3: Mean Bureaucratic Fragmentation
frag <- aggregate(frag ~ year, data = all1, FUN = function(x) mean(x, na.rm = T))
#par(mar=c(0,0,0,0))
plot(frag, type = "l", ylab = "Fragmentation (1 - HHI)", xlab = "Year", main = "Mean Bureaucratic Fragmentation")

## Figure 4: Mean Budget Total for Largest Aid Agency Types, Plurality versus PR Institutions
types <- c("Science/Education", "Economic", "Foreign", "Defense", "Development", "Domestic", "Energy")
tottype <- aggregate(ay.total ~ type, data = all1ay, FUN = sum)
tottype.pl <- aggregate(ay.total ~ type, data = subset(all1ay, pluralty == 1), FUN = mean)
tottype.pr <- aggregate(ay.total ~ type, data = subset(all1ay, pluralty == 0), FUN = mean)
tottype.pl$Plurality <- tottype.pl$ay.total
tottype.pr$PR <- tottype.pr$ay.total
tottype2 <- subset(tottype, ay.total > 5e+09)
totboth <- merge(tottype.pl[,c(1,3)], tottype.pr[,c(1,3)], by = "type")
tottype2.pl <- subset(tottype.pl, ay.total > summary(tottype.pl$ay.total)[2])
tottype2.pr <- subset(tottype.pr, ay.total > summary(tottype.pr$ay.total)[2])
tottype2.pl$Plurality <- tottype2.pl$ay.total
tottype2.pr$PR <- tottype2.pr$ay.total
totboth2 <- merge(tottype2.pl[,c(1,3)], tottype2.pr[,c(1,3)], by = "type")
totboth2$type <- ifelse(totboth2$type == "Science/Education", "Education", paste(totboth2$type))
totboth3 <- totboth2[which(totboth2$type != "Misc" & totboth2$type != "Regional"),]
par(mar=c(4,4,3,3))
barplot(t(as.matrix(totboth3[,2:3])), border = "white", names.arg = t(as.matrix(totboth3[,1])), col = c("skyblue4", "orange2"), space = 0.04, font.axis = 2, xlab = "Agency Type", main = "Mean Budget Total for Largest Aid Agency Types, \nPlurality versus PR Institutions", ylab = "Budget in Real USD")
legend(title = "Institutional Type:", "topright", colnames(totboth3[,2:3]), pt.cex = 2.5, pch = 15, col = c("skyblue4", "orange2"))

## Figure 5: Mean Bureaucratic Fragmentation in PR vs. Plurality Systems
pr <- aggregate(frag ~ year, data = subset(all1, pr == 1), FUN = mean)
pl <- aggregate(frag ~ year, data = subset(all1, pluralty == 1), FUN = mean)
prpl <- merge(pr, pl, by = "year", suffixes = c(".pr", ".pl"))
ts1 <- ts(prpl[,2:3], start = 1975)
#PLOT
par(mar=c(4,4,3,2))
plot(ts1, plot.type = "single", lty = c(2,1), col = c("orange2", "skyblue4"), xlab = "Year", ylab= "Bureaucratic Fragmentation", main = "Mean Bureaucratic Fragmentation \nin PR vs. Plurality Systems", lwd = 1.6)
legend("topleft", c("PR","Plurality"), pch = NA, lty = c(2,1), col = c("orange2", "skyblue4"), lwd = 2, title = "Institutional Type:")

## Figure 6: Mean Proportion of Aid Tied, by Agency Type
# Primary agency
tiedmain <- aggregate(tied_amount_usd_nominal ~ ay, data = subset(all2ay, main == 1), FUN = sum)
totmain <- aggregate(untied_amount_usd_nominal ~ ay, data = subset(all2ay, main == 1), FUN = sum)
percentmain <- tiedmain$tied_amount_usd_nominal/(tiedmain$tied_amount_usd_nominal+totmain$untied_amount_usd_nominal)
# Agricultural
tiedag <- aggregate(tied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Agriculture"), FUN = sum)
totag <- aggregate(untied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Agriculture"), FUN = sum)
percentag <- tiedag$tied_amount_usd_nominal/(tiedag$tied_amount_usd_nominal+totag$untied_amount_usd_nominal)
# Economic
tiedecon <- aggregate(tied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Economic"), FUN = sum)
totecon <- aggregate(untied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Economic"), FUN = sum)
percentecon <- tiedecon$tied_amount_usd_nominal/(totecon$untied_amount_usd_nominal+tiedecon$tied_amount_usd_nominal)
# Foreign
tiedfor <- aggregate(tied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Foreign" & main == 0), FUN = sum)
totfor <- aggregate(untied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Foreign" & main == 0), FUN = sum)
percentfor <- tiedfor$tied_amount_usd_nominal/(tiedfor$tied_amount_usd_nominal+totfor$untied_amount_usd_nominal)
# Science
tiedsci <- aggregate(tied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Science/Education"), FUN = sum)
totsci <- aggregate(untied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Science/Education"), FUN = sum)
percentsci <- tiedsci$tied_amount_usd_nominal/(tiedsci$tied_amount_usd_nominal+totsci$untied_amount_usd_nominal)
# Domestic
tieddom <- aggregate(tied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Domestic"), FUN = sum)
totdom <- aggregate(untied_amount_usd_nominal ~ ay, data = subset(all2ay, type == "Domestic"), FUN = sum)
percentdom <- tieddom$tied_amount_usd_nominal/(tieddom$tied_amount_usd_nominal+totdom$untied_amount_usd_nominal)
# Data for Barplot 
data <- data.frame(
  coef <- c(mean(all2ay$percent2, na.rm = T), mean(percentmain, na.rm = T), mean(percentdom, na.rm = T), mean(percentag, na.rm = T), mean(percentecon, na.rm = T), mean(percentsci, na.rm = T)),
  med <- c(median(all2ay$percent2, na.rm = T), median(percentmain, na.rm = T), median(percentdom, na.rm = T), median(percentag, na.rm = T), median(percentecon, na.rm = T), median(percentsci, na.rm = T)),
  se <- c(sd(all2ay$percent2, na.rm = T)/sqrt(length(all2ay$percent2)), sd(percentmain, na.rm = T)/sqrt(length(percentmain)), sd(percentdom, na.rm = T)/sqrt(length(percentdom)), sd(percentag, na.rm = T)/sqrt(length(percentag)), sd(percentecon, na.rm = T)/sqrt(length(percentecon)), sd(percentsci, na.rm = T)/sqrt(length(percentsci))),
  iqr <- c(IQR(all2ay$percent2, na.rm = T), IQR(percentmain, na.rm = T), IQR(percentdom, na.rm = T), IQR(percentag, na.rm = T), IQR(percentecon, na.rm = T), IQR(percentsci, na.rm = T)),
  names.arg <- c("Agencies (All)", "Primary Aid Agency", "Domestic Issues", "Agiculture", "Economic/Finance", "Education/Science")
)
# Barplot
install.packages("ggplot2")
library(ggplot2)
par(mar=c(2,3,1,1))
ggplot(data) +
  geom_bar(aes(x=names.arg, y=coef), stat="identity", fill="skyblue4", alpha=0.7) +
  geom_errorbar( aes(x=names.arg, ymin=coef-se, ymax=coef+se), width=0.4, colour="orange2", alpha=0.9, size=1.3) +
  coord_flip() +
  ggtitle("Mean Proportion of Aid Tied, by Agency Type") +
  ylab("Proportion Aid Tied") +
  xlab("") +
  theme(
    axis.text.y = element_text(size = 13),
    title = element_text(size = 15))

### APPENDIX FIGURES AND TABLES ###
### TABLES
## Table 3: Full Model of Table 1
stargazer(f1cp, f2cp, f3cp, f4cp, f5cp, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Plurality", "Veto Players", "Lib. Econ.", "Gov. Vote-Share", "Opp. Vote-Share", "HHI Government", "Fractionalization", "Leg. Election Year", "Exec. Election Year", "Unified Gov.", "Polarization", "Conservative Gov.", "GDP (trillion USD)", "GDP/cap (100,000 USD)", "Aid Budget (100 billion USD)"), 
          dep.var.labels = "Fragmentation (1 - HHI)",
          add.lines = list(c("Year FE", "N", "N", "Y", "N", "Y"),
                           c("Donor FE", "N", "N", "N", "Y", "Y"),
                           c("Observations", paste(length(f1p$fitted.values)), paste(length(f2p$fitted.values)), paste(length(f3p$fitted.values)), paste(length(f4p$fitted.values)), paste(length(f5p$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(f1p)[9]), digits = 3)), paste(round(as.numeric(summary(f2p)[9]), digits = 3)), paste(round(as.numeric(summary(f3p)[9]), digits = 3)), paste(round(as.numeric(summary(f4p)[9]), digits = 3)), paste(round(as.numeric(summary(f5p)[9]), digits = 3)))))

## Table 4: Full Model of Table 2
stargazer(tie1ayc, tie2ayc, tie4ayc, tie5ayc, tie6ayc, no.space = T, digits = 2, omit = c("year", "donor", "dy", "type"), 
          covariate.labels = c("Primary Agency", "Plurality", "Veto Players", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          add.lines = list(c("Year FE", "N", "N", "Y", "N", "Y"),
                           c("Donor FE", "N", "N", "N", "Y", "Y"),
                           c("Observations", paste(length(tt1$fitted.values)), paste(length(tt2$fitted.values)), paste(length(tt4$fitted.values)), paste(length(tt5$fitted.values)), paste(length(tt6$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(tt1)[9]), digits = 3)), paste(round(as.numeric(summary(tt2)[9]), digits = 3)), paste(round(as.numeric(summary(tt4)[9]), digits = 3)), paste(round(as.numeric(summary(tt5)[9]), digits = 3)), paste(round(as.numeric(summary(tt6)[9]), digits = 3)))))

## Table 5: Personal-Vote Models
f1 <- lm(frag ~ pers.rank, data = all1)
f2 <- lm(frag ~ pers.rank + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f3 <- lm(frag ~ factor(year) + pers.rank + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f4 <- lm(frag ~ factor(year) + countryname + pers.rank + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
vcov1f <- vcovCL(f1, cluster = all1$year)
vcov2f <- vcovCL(f2, cluster = all1$year)
vcov3f <- vcovCL(f3, cluster = all1$year)
vcov4f <- vcovCL(f4, cluster = all1$year)
f1c <- coeftest(f1, vcov = vcov1f)
f2c <- coeftest(f2, vcov = vcov2f)
f3c <- coeftest(f3, vcov = vcov3f)
f4c <- coeftest(f4, vcov = vcov4f)
stargazer(f1c, f2c, f3c, f4c, star.cutoffs = c(0.05, 0.01), no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Personal Vote", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          add.lines = list(c("Year FE", "N", "N", "Y", "Y", "N", "N", "Y", "Y"),
                           c("Donor FE", "N", "N", "N", "Y", "N", "N", "N", "Y"),
                           c("Observations", paste(length(f1$fitted.values)), paste(length(f2$fitted.values)), paste(length(f3$fitted.values)), paste(length(f4$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(f1)[9]), digits = 3)), paste(round(as.numeric(summary(f2)[9]), digits = 3)), paste(round(as.numeric(summary(f3)[9]), digits = 3)), paste(round(as.numeric(summary(f4)[9]), digits = 3)))))

## Table 6: Number of Aid Agencies as DV
summary(f1n <- lm(number ~ pluralty, data = all1))
summary(f2n <- lm(number ~ pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
summary(f3n <- lm(number ~ factor(year) + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
summary(f4n <- lm(number ~ factor(year) + countryname + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
vcov1n <- vcovCL(f1n)
vcov2n <- vcovCL(f2n, cluster = all1$countryname)
vcov3n <- vcovCL(f3n, cluster = all1$countryname)
vcov4n <- vcovCL(f4n, cluster = all1$countryname)
t1n <- coeftest(f1n, vcov = vcov1n)
t2n <- coeftest(f2n, vcov = vcov2n)
t3n <- coeftest(f3n, vcov = vcov3n)
t4n <- coeftest(f4n, vcov = vcov4n)
stargazer(t1n, t2n, t3n, t4n, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Plurality", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          dep.var.labels = "Number of Aid Agencies",
          add.lines = list(c("Year FE", "N", "N", "Y", "Y"),
                           c("Donor FE", "N", "N", "N", "Y"),
                           c("Observations", paste(length(f1n$fitted.values)), paste(length(f2n$fitted.values)), paste(length(f3n$fitted.values)), paste(length(f4n$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(f1n)[9]), digits = 3)), paste(round(as.numeric(summary(f2n)[9]), digits = 3)), paste(round(as.numeric(summary(f3n)[9]), digits = 3)), paste(round(as.numeric(summary(f4n)[9]), digits = 3)))))

## Table 7: Recipient-Level Fragmentation as DV
f1r <- lm(frag.recip ~ pluralty, data = all1)
f2r <- lm(frag.recip ~ pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f3r <- lm(frag.recip ~ factor(year) + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f4r <- lm(frag.recip ~ factor(year) + countryname + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
vcov1fr <- vcovCL(f1r, cluster = all1$countryname)
vcov2fr <- vcovCL(f2r, cluster = all1$countryname)
vcov3fr <- vcovCL(f3r, cluster = all1$countryname)
vcov4fr <- vcovCL(f4r, cluster = all1$countryname)
f1cr <- coeftest(f1r, vcov = vcov1fr)
f2cr <- coeftest(f2r, vcov = vcov2fr)
f3cr <- coeftest(f3r, vcov = vcov3fr)
f4cr <- coeftest(f4r, vcov = vcov4fr)
stargazer(f1cr, f2cr, f3cr, f4cr, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Plurality Rule", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          add.lines = list(c("Year FE", "N", "N", "Y", "Y"),
                           c("Donor FE", "N", "N", "N", "Y"),
                           c("Observations", paste(length(f1r$fitted.values)), paste(length(f2r$fitted.values)), paste(length(f3r$fitted.values)), paste(length(f4r$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(f1r)[9]), digits = 3)), paste(round(as.numeric(summary(f2r)[9]), digits = 3)), paste(round(as.numeric(summary(f3r)[9]), digits = 3)), paste(round(as.numeric(summary(f4r)[9]), digits = 3)))))

## Table 8: Limited Donor-Year Sample
summary(f1a <- lm(frag ~ pluralty, data = all))
summary(f2a <- lm(frag ~ pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all))
summary(f3a <- lm(frag ~ factor(year) + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all))
summary(f4a <- lm(frag ~ factor(year) + countryname + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all))
summary(f5a <- lm(frag ~ factor(year) + countryname + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all))
vcov1fa <- vcovCL(f1a, cluster = all$countryname)
vcov2fa <- vcovCL(f2a, cluster = all$countryname)
vcov3fa <- vcovCL(f3a, cluster = all$countryname)
vcov4fa <- vcovCL(f4a, cluster = all$countryname)
f1ca <- coeftest(f1a, vcov = vcov1fa)
f2ca <- coeftest(f2a, vcov = vcov2fa)
f3ca <- coeftest(f3a, vcov = vcov3fa)
f4ca <- coeftest(f4a, vcov = vcov4fa)
stargazer(f1ca, f2ca, f3ca, f4ca, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Plurality Rule", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          add.lines = list(c("Year FE", "N", "N", "Y", "Y"),
                           c("Donor FE", "N", "N", "N", "Y"),
                           c("Observations", paste(length(f1a$fitted.values)), paste(length(f2a$fitted.values)), paste(length(f3a$fitted.values)), paste(length(f4a$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(f1a)[9]), digits = 3)), paste(round(as.numeric(summary(f2a)[9]), digits = 3)), paste(round(as.numeric(summary(f3a)[9]), digits = 3)), paste(round(as.numeric(summary(f4a)[9]), digits = 3)))))

## Table 9: Limited Agency-Year Sample
tt1a <- lm(percent2 ~ main, data = aylimit)
tt2a <- lm(percent2 ~ main + pers.rank + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = aylimit)
tt4a <- lm(percent2 ~ factor(year) + main + pers.rank + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = aylimit)
tt5a <- lm(percent2 ~ factor(year) + donor + main + pers.rank + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = aylimit)
vcov1a <- vcovCL(tt1a, cluster = aylimit$countryname)
vcov2a <- vcovCL(tt2a, cluster = aylimit$countryname)
vcov4a <- vcovCL(tt4a, cluster = aylimit$countryname)
vcov5a <- vcovCL(tt5a, cluster = aylimit$countryname)
tie1ayca <- coeftest(tt1a, vcov = vcov1a)
tie2ayca <- coeftest(tt2a, vcov = vcov2a)
tie4ayca <- coeftest(tt4a, vcov = vcov4a)
tie5ayca <- coeftest(tt5a, vcov = vcov5a)
stargazer(tie1ayca, tie2ayca, tie4ayca, tie5ayca, no.space = T, digits = 2, omit = c("year", "donor", "dy", "type"), 
          covariate.labels = c("Primary Agency", "Personal Vote", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          add.lines = list(c("Year FE", "N", "N", "Y", "Y"),
                           c("Donor FE", "N", "N", "Y", "N"),
                           c("Agency FE", "N", "N", "N", "Y"),
                           c("Observations", paste(length(tt1a$fitted.values)), paste(length(tt2a$fitted.values)), paste(length(tt4a$fitted.values)), paste(length(tt5a$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(tt1a)[9]), digits = 3)), paste(round(as.numeric(summary(tt2a)[9]), digits = 3)), paste(round(as.numeric(summary(tt4a)[9]), digits = 3)), paste(round(as.numeric(summary(tt5a)[9]), digits = 3)))))

## Table 10: Interacting Plurality and Election Year
f1el <- lm(frag ~ exelec, data = all1)
f1eli <- lm(frag ~ pluralty*exelec, data = all1)
f2eli <- lm(frag ~ pluralty*exelec + lme + gov1vote + opp1vote + herfgov + frac + legelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f3eli <- lm(frag ~ factor(year) + pluralty*exelec + lme + gov1vote + opp1vote + herfgov + frac + legelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
f4eli <- lm(frag ~ factor(year) + countryname + pluralty*exelec + lme + gov1vote + opp1vote + herfgov + frac + legelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1)
vcov1el <- vcovCL(f1el, cluster = all1$countryname)
vcov1eli <- vcovCL(f1eli, cluster = all1$countryname)
vcov2eli <- vcovCL(f2eli, cluster = all1$countryname)
vcov3eli <- vcovCL(f3eli, cluster = all1$countryname)
vcov4eli <- vcovCL(f4eli, cluster = all1$countryname)
f1cel <- coeftest(f1el, vcov = vcov1el)
f1celi <- coeftest(f1eli, vcov = vcov1eli)
f2celi <- coeftest(f2eli, vcov = vcov2eli)
f3celi <- coeftest(f3eli, vcov = vcov3eli)
f4celi <- coeftest(f4eli, vcov = vcov4eli)
stargazer(f1cel, f1celi, f2celi, f3celi, f4celi, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Plurality", "Exec. Election Year", "Lib. Econ.", "Gov. Vote-Share", "Opp. Vote-Share", "HHI Government", "Fractionalization", "Leg. Election Year", "Unified Gov.", "Polarization", "Conservative Gov.", "GDP (trillion USD)", "GDP/cap (100,000 USD)", "Aid Budget (100 billion USD)"), 
          dep.var.labels = "Fragmentation (1 - HHI)",
          add.lines = list(c("Year FE", "N", "N", "Y", "N"),
                           c("Donor FE", "N", "N", "N", "Y"),
                           c("Observations", paste(length(f1el$fitted.values)), paste(length(f1eli$fitted.values)), paste(length(f2eli$fitted.values)), paste(length(f3eli$fitted.values)), paste(length(f4eli$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(f1el)[9]), digits = 3)), paste(round(as.numeric(summary(f1eli)[9]), digits = 3)), paste(round(as.numeric(summary(f2eli)[9]), digits = 3)), paste(round(as.numeric(summary(f3eli)[9]), digits = 3)), paste(round(as.numeric(summary(f4eli)[9]), digits = 3)))))

## Table 11: Affective Polarization
### Polarized public leads to less fragmentation (because powerful parties)
summary(polar1 <- lm(frag ~ partisan_affect_polarization, data = all1))
summary(polar2 <- lm(frag ~ partisan_affect_polarization + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
summary(polar4 <- lm(frag ~ partisan_affect_polarization + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2 + factor(year), data = all1))
summary(polar5 <- lm(frag ~ partisan_affect_polarization + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2 + countryname + factor(year), data = all1))
vcov1po <- vcovCL(polar1, cluster = all1$countryname)
vcov2po <- vcovCL(polar2, cluster = all1$countryname)
vcov4po <- vcovCL(polar4, cluster = all1$countryname)
vcov5po <- vcovCL(polar5, cluster = all1$countryname)
f1po <- coeftest(polar1, vcov = vcov1po)
f2po <- coeftest(polar2, vcov = vcov2po)
f4po <- coeftest(polar4, vcov = vcov4po)
f5po <- coeftest(polar5, vcov = vcov5po)
stargazer(f1po, f2po, f4po, f5po, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Polarization", "Plurality Rule", "Lib. Econ.", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget"), 
          add.lines = list(c("Year FE", "N", "N", "Y", "Y"),
                           c("Donor FE", "N", "N", "N", "Y"),
                           c("Observations", paste(length(polar1$fitted.values)), paste(length(polar2$fitted.values)), paste(length(polar4$fitted.values)), paste(length(polar5$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(polar1)[9]), digits = 3)), paste(round(as.numeric(summary(polar2)[9]), digits = 3)), paste(round(as.numeric(summary(polar4)[9]), digits = 3)), paste(round(as.numeric(summary(polar5)[9]), digits = 3)))))

## Table 12: Interacting Bureaucracy with Plurality
pl <- subset(all2ay, pluralty == 1)
pr <- subset(all2ay, pluralty == 0)
ttp <- lm(pct.rescale ~ factor(year) + donor + main*number + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = pl)
ttnp <- lm(pct.rescale ~ factor(year) + donor + main*number + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = pr)
ttp <- lm(pct.rescale ~ factor(year) + donor + main*number + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = pl)
ttnp <- lm(pct.rescale ~ factor(year) + donor + main*number + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = pr)
tt1 <- lm(pct.rescale ~ main*number, data = all2ay)
tt5 <- lm(pct.rescale ~ factor(year) + donor + main*number + pluralty + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all2ay)
vcovp <- vcovCL(ttp, cluster = pl$donor)
vcovnp <- vcovCL(ttnp, cluster = pr$donor)
vcov1 <- vcovCL(tt1, cluster = all2ay$donor)
vcov5 <- vcovCL(tt5, cluster = all2ay$donor)
tie1ayc <- coeftest(tt1, vcov = vcov1)
tie5ayc <- coeftest(tt5, vcov = vcov5)
tiep <- coeftest(ttp, vcov = vcovp)
tienp <- coeftest(ttnp, vcov = vcovnp)
stargazer(tie1ayc, tie5ayc, tiep, tienp, no.space = T, digits = 2, omit = c("year", "donor", "dy", "type"), 
          covariate.labels = c("Primary Agency", "Num. of Agencies", "Plurality", "Gov. Vote", "Opp. Vote", "HHI Gov", "Frac.", "Leg Elec", "Exec Elec", "Unified Gov.", "Polariz,", "Cons.", "GDP", "GDP/cap", "Aid Budget", "PrimaryxNumber"), 
          add.lines = list(c("Year FE", "N", "Y", "Y", "Y"),
                           c("Donor FE", "N", "Y", "Y", "Y"),
                           c("Observations", paste(length(tt1$fitted.values)), paste(length(tt5$fitted.values)), paste(length(ttp$fitted.values)), paste(length(ttnp$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(tt1)[9]), digits = 3)), paste(round(as.numeric(summary(tt5)[9]), digits = 3)), paste(round(as.numeric(summary(ttp)[9]), digits = 3)), paste(round(as.numeric(summary(ttnp)[9]), digits = 3)))))

## Table 13: Subsets of Institutional Variables
summary(sub2 <- lm(frag ~ pluralty + checks + lme + pers.rank + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
summary(sub3 <- lm(frag ~ pluralty + checks + pers.rank + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
summary(sub4 <- lm(frag ~ pluralty + lme + pers.rank + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
summary(sub5 <- lm(frag ~ pluralty + checks + lme + gov1vote + opp1vote + herfgov + frac + legelec + exelec + allhouse + polariz + conservatism + gdp2 + pcap2 + total2, data = all1))
vcov2sub <- vcovCL(sub2, cluster = all1$countryname)
vcov3sub <- vcovCL(sub3, cluster = all1$countryname)
vcov4sub <- vcovCL(sub4, cluster = all1$countryname)
vcov5sub <- vcovCL(sub5, cluster = all1$countryname)
sub2p <- coeftest(sub2, vcov = vcov2sub)
sub3p <- coeftest(sub3, vcov = vcov3sub)
sub4p <- coeftest(sub4, vcov = vcov4sub)
sub5p <- coeftest(sub5, vcov = vcov5sub)
stargazer(sub2p, sub3p, sub4p, sub5p, no.space = T, digits = 2, omit = c("year", "countryname"), 
          covariate.labels = c("Plurality", "Veto Players", "Lib. Econ.", "Personal Vote", "Gov. Vote-Share", "Opp. Vote-Share", "HHI Government", "Fractionalization", "Leg. Election Year", "Exec. Election Year", "Unified Gov.", "Polarization", "Conservative Gov.", "GDP (trillion USD)", "GDP/cap (100,000 USD)", "Aid Budget (100 billion USD)"), 
          dep.var.labels = "Fragmentation (1 - HHI)",
          add.lines = list(c("Observations", paste(length(sub2$fitted.values)), paste(length(sub3$fitted.values)), paste(length(sub4$fitted.values)), paste(length(sub5$fitted.values))),
                           c("Mult. R-2", paste(round(as.numeric(summary(sub2)[9]), digits = 3)), paste(round(as.numeric(summary(sub3)[9]), digits = 3)), paste(round(as.numeric(summary(sub4)[9]), digits = 3)), paste(round(as.numeric(summary(sub5)[9]), digits = 3)))))

## Table 14: Summary Statistics
stargazer(all1[,c("frag", "pers.rank", "lme", "gov1vote", "opp1vote", "herfgov", "frac", "legelec", "exelec", "allhouse", "polariz", "conservatism", "gdp2", "pcap2", "total2")], covariate.labels = c("Fragmentation", "Personal Vote Index", "Lib. Mark. Econ.", "Gov. Vote-Share", "Opp. Vote-Share", "HHI Gov", "Frac.", "Elec. Year (Leg)", "Elec. Year (Exec)", "Unified Gov.", "Polarization", "Conservativism", "GDP (bil USD)", "GDP/cap (100,000)", "Aid Budget (bil USD)"))
stargazer(all2ay[,c("percent2", "pers.rank", "lme", "gov1vote", "opp1vote", "herfgov", "frac", "legelec", "exelec", "allhouse", "polariz", "conservatism", "gdp2", "pcap2", "total2")], covariate.labels = c("Percent Tied Aid", "Personal Vote Index", "Lib. Mark. Econ.", "Gov. Vote-Share", "Opp. Vote-Share", "HHI Gov", "Frac.", "Elec. Year (Leg)", "Elec. Year (Exec)", "Unified Gov.", "Polarization", "Conservativism", "GDP (bil USD)", "GDP/cap (100,000)", "Aid Budget (bil USD)"))

## Table 15: Correlation Matrices
stargazer(cor(f2p$model))
stargazer(cor(tt2$model))

### APPENDIX FIGURES
## Figure 7: Correlation of bureaucratic and recipient-level fragmentation among major aid donors
install.packages("calibrate")
library(calibrate)
frag <- aggregate(all1$frag ~ all1$countryname, FUN = mean)
frag.recip <- aggregate(all1$frag.recip ~ all1$countryname, FUN = mean)
allfrag <- merge(frag, frag.recip, by = "all1$countryname")
plot(allfrag$`all1$frag`, allfrag$`all1$frag.recip`, xlab = "Organizational Fragmentation", ylab = "Fragmentation among Recipients", main = "Measures of Aid Fragmentation", )
textxy(allfrag$`all1$frag`, allfrag$`all1$frag.recip`, allfrag$`all1$countryname`, cex = .7, pos = 4)

## Figure 8: Marginal effects plots of number of agencies and aid tying
install.packages("interplot")
library(interplot)
ttp <- lm(pct.rescale ~ main*number, data = subset(all2ay, pluralty == 1))
ttnp <- lm(pct.rescale ~ main*number, data = subset(all2ay, pluralty == 0))
interplot(ttp, "main", "number", hist = T) +
  xlab("Number of Agencies") +
  ylab("Marginal Effects of Lead Agency") +
  ggtitle("Marginal Effects Plot, Plurality Institutions")
interplot(ttnp, "main", "number", hist = T) +
  xlab("Number of Agencies") +
  ylab("Marginal Effects of Lead Agency") +
  ggtitle("Marginal Effects Plot, PR Institutions")

